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Abstract 

The  recently  developed  time-periodic  fluctuation-dissipation  theorem  (FDT) 
provides  a  very  convenient  way  of  addressing  the  climate  change  of  atmospheric 
systems  with  seasonal  cycle  by  utilizing  statistics  of  the  present  climate.  A 
triad  nonlinear  stochastic  model  with  exactly  solvable  first  and  second  order 
statistics  is  introduced  here  as  an  unambiguous  test  model  for  FDT  in  a  time- 
periodic  setting.  This  model  mimics  the  nonlinear  interaction  of  two  Rossby 
waves  forced  by  baroclinic  processes  with  a  zonal  jet  forced  by  a  polar  temper¬ 
ature  gradient.  Periodic  forcing  naturally  introduces  the  seasonal  cycle  into  the 
model.  The  exactly  solvable  first  and  second  order  statistics  are  utilized  to  com¬ 
pute  both  the  ideal  mean  and  variance  response  to  the  perturbations  in  forcing 
or  dissipation  and  the  quasi-Gaussian  approximation  of  FDT  (qG-FDT)  that 
uses  the  mean  and  the  covariance  in  the  equilibrium  state.  The  time-averaged 
mean  and  variance  qG-FDT  response  to  perturbations  of  forcing  or  dissipation 
is  compared  with  the  corresponding  ideal  response  utilizing  the  triad  test-model 
in  a  number  of  regimes  with  various  dynamical  and  statistical  properties  such 
as  weak  or  strong  non-Gaussianity  and  resonant  or  non-resonant  forcing.  It  is 
shown  that  even  in  a  strongly  non-Gaussian  regime,  qG-FDT  has  surprisingly 
high  skill  for  the  mean  response  to  the  changes  in  forcing.  On  the  other  hand, 
the  performance  of  qG-FDT  for  the  variance  response  to  the  perturbations  of 
dissipation  is  good  in  the  near-Gaussian  regime  and  deteriorates  in  the  strongly 
non-Gaussian  regime.  The  results  here  on  the  test  model  should  provide  useful 
guidelines  for  applying  the  time-periodic  FDT  to  more  complex  realistic  systems 
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1.  Introduction 

The  application  of  the  fluctuation-dissipation  theorem  (FDT)  [5,  28]  to  cli¬ 
mate  change  science  has  a  growing  interest  among  researchers.  The  classical 
FDT  roughly  states  that  in  order  to  obtain  the  mean  response  of  the  system 
of  identical  particles  in  statistical  equilibrium  to  small  external  perturbations, 
it  is  sufficient  to  find  certain  correlation  functions  of  the  unperturbed  system. 
Motivated  by  Kraichnan’s  generalization  of  FDT  to  systems  with  the  Liouville 
property  [19,  20],  Leith  [21]  suggested  that  if  the  climate  system  satisfied  a 
suitable  fluctuation-dissipation  theorem  (FDT)  then  climate  response  to  small 
external  forcing  or  other  parameter  changes  could  be  calculated  by  estimating 
suitable  statistics  in  the  present  climate.  For  the  general  FDT,  see  [9],  [31],  and 
[23].  The  important  practical  and  conceptual  advantages  for  climate  change 
science  when  a  skillful  FDT  algorithm  can  be  established  is  that  the  linear  sta¬ 
tistical  response  operator  produced  by  FDT  can  be  utilized  directly  for  multiple 
climate  change  scenarios,  multiple  changes  in  forcing,  dissipation  and  other  pa¬ 
rameters,  and  inverse  modelling  directly  [14,  16]  without  the  need  of  running  the 
complex  climate  model  in  each  individual  case,  often  a  computational  problem  of 
overwhelming  complexity.  With  these  interesting  possibilities  for  FDT,  Leith’s 
suggestion  inspired  a  first  wave  of  systematic  research  [6,  7,  29,  17,  13,  15,  14,  16] 
for  various  idealized  climate  models  for  the  mean  response  to  changes  in  exter¬ 
nal  forcing.  All  of  this  work  utilizes  the  quasi-Gaussian  approximation  (qG- 
FDT)  suggested  by  Leith  [21,  23,  4].  Recently,  mathematical  theory  for  FDT 
[23,  27,  25,  22]  has  supplied  important  generalizations  and  new  ways  to  interpret 
FDT.  These  developments  have  lead  to  improved  theoretical  understanding  of 
the  qG-FDT  algorithm,  new  applications  of  qG-FDT  beyond  the  mean  response 
with  significant  skill  [16]  and  new  computational  algorithms  beyond  qG-FDT 
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with  improved  high  skill  for  both  the  mean  and  variance  of  low  frequency  re¬ 
sponse  which  have  been  tested  in  a  variety  of  models  [2,  3,  4] .  There  is  also  recent 
mathematical  theory  for  identifying  the  “most  dangerous”  perturbations  in  a 
given  class  [23,  27].  The  above  work  on  algorithms  for  FDT  in  climate  change 
science  is  for  forced  dissipative  large  dimensional  dynamical  systems.  Thus,  the 
attractor  is  a  fractal  object,  typically,  but  the  presence  of  large  dimensional 
unstable  manifolds  helps  FDT  algorithms  (see  [28,  8])  to  overcome  this  lack  of 
smoothness  for  these  high  dimensional  systems.  Furthermore,  the  short-time 
FDT  algorithms  developed  in  [2,  3,  4]  work  for  systems  without  any  smooth¬ 
ness  on  the  attractor.  However,  most  of  the  theories  and  algorithms  developed 
for  FDT  assume  a  stationary  climate,  which  excludes  the  study  of  important 
practical  issues  of  climate  change  science  involving  time-periodic  statistics  such 
as  the  diurnal  or  seasonal  cycle. 

Very  recently  generalizations  of  FDT  to  time-dependent  ensembles  were  de¬ 
veloped  by  Majda  and  Wang  [27].  As  a  special  case  of  more  general  results 
in  their  paper,  Majda  and  Wang  developed  the  FDT  for  time-periodic  systems 
as  well  as  approximate  algorithms  for  climate  response  with  a  seasonal  cycle. 
In  particular,  a  quasi-Gaussian  algorithm  for  computing  FDT  response  of  a 
stochastic  system  with  time-periodic  invariant  measure  was  proposed  in  [27].  In 
a  time-periodic  setting,  the  FDT  helps  to  answer  such  practical  questions  as  how 
will  the  time  averaged  monthly,  seasonal,  or  yearly  mean  or  variance  of  certain 
physical  variables  change  if  the  external  forcing  is  perturbed.  For  example,  one 
can  be  interested  in  how  much  the  averaged  temperature  in  the  month  of  April 
will  change  if  the  forcing  becomes  stronger  in  the  month  of  January. 

In  this  paper,  we  propose  a  simple  yet  very  rich  triad  test  model  with  both 
periodic  deterministic  forcing  and  stochastic  forcing  as  a  test  bed  for  time  pe¬ 
riodic  FDT.  This  model  has  already  been  introduced  and  used  by  the  authors 
in  applications  to  filtering  problems  with  multiple  time  scales  [11,  12].  The 
presence  of  noise  makes  the  attractor  smooth  and  there  is  a  general  rigorous 
justification  of  FDT  in  this  context  [18].  The  mathematical  formulation  of  FDT 
as  linear  response  theory  for  forced  dissipative  stochastic  dynamical  systems  is 
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an  appropriate  setting  for  these  applications  since  many  current  improvements 
in  the  comprehensive  computer  models  for  climate  change  involve  stochastic 
components  [30],  while  lower  dimensional  reduced  models  typically  also  involve 
stochastic  noise  terms  [24,  22,  25]. 

The  key  properties  of  the  test  model  which  make  it  very  attractive  for  testing 
FDT  are 

•  its  exactly  solvable  structure,  i.e. ,  mean,  variance  and  in  principle  any 
higher  order  moments  can  be  computed  analytically,  therefore,  the  full 
ideal  response  to  external  perturbations  can  be  found  exactly;  further¬ 
more,  the  qG-FDT  algorithm  can  be  applied  and  tested  because  it  re¬ 
quires  knowledge  of  the  time-periodic  mean  and  covariance  statistics  of 
the  unperturbed  climate, 

•  natural  time-periodic  forcing  to  study  the  performance  of  FDT  on  the 
systems  with  seasonal  cycle, 

•  nonlinear  dynamics  which  provides  an  opportunity  to  study  both  Gaus¬ 
sian  and  non-Gaussian  regimes  with  the  quasi-Gaussian  approximation  of 
FDT. 

The  model  is  a  triad  nonlinear  stochastic  model  consisting  of  one  real  mode,  u\, 
and  one  complex  mode,  it 2,  that  interacts  with  u\  through  catalytic  nonlinear 
coupling 

=  -71w1  +  fi(t)  +  aiWu  (1) 

at 

—jt  =  (-72  +i(wo  +  aoUi))u2  +  /2W  +  02W2,  (2) 

at 

where  fi{t)  and  /2(f)  are  periodic  functions  of  time  with  the  same  period  To, 
which  represents  the  annual  cycle,  71,  72,  07,  and  (J2  are  dissipation  and  stochas¬ 
tic  forcing  coefficients  that  represent  the  interaction  of  the  triad  system  with 
other  unresolved  modes,  ujq  is  deterministic  part  of  the  frequency  for  the  com¬ 
plex  mode,  u 2,  and  a 0  is  the  coefficient  measuring  nonlinearity.  In  [11,  12]  the 
authors  demonstrated  how  to  use  the  special  structure  of  the  nonlinearity  to  find 
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analytical  formulas  for  the  first  and  second  order  statistics  of  system  (1),  (2). 
Here,  we  use  the  long  time  limit  of  those  analytical  formulas  in  order  to  find 
statistics  in  a  time-periodic  statistical  equilibrium.  Then,  these  first  and  second 
order  statistics  can  be  used  for  computing 

•  the  “ideal”  mean  and  variance  response  of  the  system  to  the  perturbations 
of  external  forcing  and  dissipation, 

•  the  quasi-Gaussian  approximation  of  FDT,  which  requires  the  knowledge 
of  mean  and  covariance  statistics  of  the  unperturbed  climate. 

The  test  model  (1),  (2)  is  motivated  by  the  interaction  of  a  barotropic  or  baro- 
clinic  Rossby  wave  on  a  sphere  represented  by  mode  the  U2  with  Rossby  fre¬ 
quency  to o  with  a  strong  zonal  wind  represented  by  mode  U\ .  A  special  solution 
to  a  model  similar  to  system  (1),  (2)  was  given  in  [26]  but  for  the  system  without 
stochastic  forcing.  The  forcing  fi(t)  represents  the  direct  forcing  of  the  zonal 
jet  from  the  polar  temperature  gradient.  On  the  other  hand,  /2(f)  models  the 
forcing  of  Rossby  waves  due  to  baroclinic  moist  processes  or  sea  surface  temper¬ 
ature.  Naturally,  both  of  these  components  of  external  forcing  have  a  seasonal 
cycle.  The  advection  of  the  Rossby  wave  by  the  zonal  jet  is  modeled  by  the 
nonlinear  term  with  coupling  coefficient  a o-  By  varying  the  parameters  of  the 
test  model  we  can  mimic  different  scenarios  with  various  ratios  of  the  energies 
of  u\  and  112,  various  characteristic  time  scales  of  these  modes  compared  to  the 
seasonal  cycle  To,  various  Rossby  wave  frequencies  when  compared  with  the 
seasonal  cycle  (external  forcing)  frequency,  and  also  various  strengths  of  nonlin¬ 
earity.  By  varying  the  nonlinearity  strength,  ao,  we  can  control  the  departure 
of  the  system  statistics  from  the  Gaussian  state.  As  in  climate  science,  we  call 
the  period  [0,  To]  a  “year”  that  consists  of  four  equal  “seasons”  and  each  season 
is  divided  into  three  equal  “months” . 

We  compute  the  the  time-averaged  mean  response  to  the  changes  in  forcing 
and  variance  response  to  the  changes  in  dissipation  using  qG-FDT  algorithms. 
Then,  we  compare  these  qG-FDT  responses  that  are  given  by  the  linear  op¬ 
erators  with  the  corresponding  ideal  response  operators  that  can  be  computed 
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using  exactly  solvable  statistics.  We  choose  a  number  of  test-cases  with  various 
types  of  system  behavior  to  test  the  skill  of  the  qG-FDT  algorithm  for  the  triad 
test  model  (35),  (36).  We  will  start  by  studying  near-Gaussian  regime  where 
the  skill  of  qG-FDT  is  very  high.  Then,  we  will  consider  an  interesting  test-case, 
when  the  statistics  averaged  over  a  certain  month  or  over  a  full  year  are  near- 
Gaussian,  while  the  averages  over  a  specific  season  are  strongly  non-Gaussian. 
In  a  near-Gaussian  regime,  we  expect  high  skill  of  the  quasi-Gaussian  approxi¬ 
mation  of  FDT  whereas  the  strongly  non-Gaussian  regime  is  a  tough  test  case 
for  quasi-Gaussian  approximation  and  its  skill  can  deteriorate.  Moreover,  we 
will  study  both  resonant  and  non-resonant  situations,  when  the  Rossby  wave 
frequency,  uq  is  either  equal  or  different  from  the  external  forcing  frequency.  We 
will  find  that  even  though  in  the  resonant  case,  the  system  becomes  strongly 
nonlinear  and  non-Gaussian,  the  skill  of  qG-FDT  for  the  mean  response  to  the 
changes  in  forcing  is  surprisingly  high  and  comparable  to  the  skill  of  qG-FDT 
in  the  non-resonant  and  more  Gaussian  regimes.  On  the  other  hand,  the  skill 
of  qG-FDT  for  the  variance  response  to  the  changes  of  dissipation  deteriorates 
significantly  as  the  nonlinearity  increases. 

The  rest  of  the  paper  is  organized  as  follows.  In  Section  2,  we  briefly  summa¬ 
rize  the  general  theory  for  time  dependent  FDT.  In  Section  3,  we  give  a  detailed 
discussion  of  the  triad  model  (1),  (2),  its  solution  and  time-periodic  equilib¬ 
rium  mean  and  covariance.  There,  we  also  demonstrate  how  to  use  the  exactly 
solvable  mean  and  covariance  in  order  to  find  the  ideal  mean  and  variance  re¬ 
sponse  to  the  changes  in  forcing  and  dissipation.  In  Section  4,  we  compute 
quasi-Gaussian  approximation  to  the  mean  and  variance  response  of  the  triad 
system  (1),  (2)  to  the  perturbations  of  forcing  or  dissipation.  In  Section  5,  we 
present  the  results  of  our  study  of  the  skill  of  the  qG-FDT  for  the  triad  system. 
Finally,  in  Section  6,  we  summarize  the  results  of  the  paper  and  discuss  future 
work. 
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2.  Theory  for  time-periodic  FDT 

In  this  Section,  we  briefly  summarize  the  theory  for  the  time-dependent  FDT 
in  a  time-periodic  statistical  steady  state.  A  much  more  extensive  discussion 
of  the  subject  can  be  found  in  [27].  Consider  a  generic  well-posed  system  of 
Stochastic  Differential  Equations  (SDE)  in  the  Ito  form 

nil 

—  =  F(u,t)  +  a(u,  t)W (t),  (3) 

at 

where  u  €  RN,  F  £  RN ,  W  is  M-dimensional  white  noise  in  time,  and  cr  is  an 
N  x  M  matrix.  We  assume  that  both  F(u,t)  and  a(u,t )  are  periodic  functions 
of  time  with  the  same  period  To,  i.e. , 

F(u,t  +  T0)  =  F(u,t ), 

a(u,t  +  T0)  =  a{u,t). 

Equation  (3)  models  the  motion  of  some  physical  system.  Because  of  the  time- 
dependent  forcing  and  noise,  system  (3)  does  not  reach  any  time-independent 
statistical  equilibrium.  However,  we  can  consider  time-dependent  statistical 
equilibrium  of  this  system.  Of  course,  even  the  time-dependent  statistical  equi¬ 
librium  may  not  exist  for  an  arbitrary  system  (3).  However,  if  the  system  is  dis¬ 
sipative  in  certain  appropriate  sense  the  existence  of  the  statistical  equilibrium 
can  be  established  ([27]  and  references  therein).  In  particular,  we  assume  that 
the  time-periodic  equilibrium  is  described  by  the  time-periodic  pdf,  peq(u,s), 
with  Peqiu ,  s  +  To)  =  Peqiu ,  s)  which  satisfies  the  Fokker-Planck  equation 

“  V“  '  '  Vu(aaTPeq)  =  0.  (4) 

Naturally,  two  types  of  averaging  arise 

•  phase  average:  for  any  function  G(u,  s),  we  have 

(G(u,s))(s)  =  J  G(u,s)peq(u,s)du,  (5) 

•  time  average:  for  any  periodic  function  f(s) 

(/(s))t0  =  -7-  /  f{s)ds.  (6) 

Jo  J o 
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The  combined  average  over  both  phase  space  and  time  is  defined  as 


~(G(u,s))  =  {{G(u,s)))To-  (7) 

Note  that  with  such  averaging,  peq(u ,  s)  becomes  a  probability  measure  on  the 
space  x  S1 

jr  J  j Peq{u,s)duds  =  1.  (8) 

Suppose,  we  are  interested  in  how  the  mean  of  some  nonlinear  functional, 
A(u,  s),  changes  when  a  small  perturbation  is  applied  to  the  forcing  F(u,  t).  We 
consider  the  perturbations  of  the  general  type 

6F(u,t)  =  a(u,s)Sf(t),  (9) 

where  a(u,  s)  is  time-periodic  vector  function,  f(t)  is  some  scalar  function  of 
time,  and  5  is  a  small  parameter.  Then  the  time-periodic  FDT  states  [27]  that 
the  finite  time  response  of  the  mean  of  the  nonlinear  functional  A(u,  s)  after 
time  t  is  given  by 

S(A(u,s))=  [  R(t  —  t')Sf(t')dt',  (10) 

Jo 

where  the  response  operator  is  computed  via 

R(t)  =  ~(A(u(t  +  s),  t  +  s)  <g>  Bp(u(s),  s)),  (11) 

and  u(t)  is  the  solution  of  the  phase-shifted  dynamical  equation 
du 

—  =  F(u,t  +  s)  +  a(u,t  +  s)W,  (12) 

dt 

with  the  initial  condition  given  as  the  value  of  the  trajectory  at  time  s 

u\t=o  =  u(s).  (13) 

The  functional  Bp{u,s)  in  (11)  has  the  explicit  form 

Vn(a(«,  S)peq(u,  s)) 

Peg(u,s) 

8 


Bf(u,s ) 


(14) 


The  infinite  time  response  of  the  system  to  the  time-independent  perturbation 
with  a  constant  change  in  forcing,  <5/(f)  =  const,  is  given  by  the  linear  relation¬ 
ship 

5(A(u,s)j  =  U6f,  (15) 


where  the  response  operator  has  the  form 


n=  R(t)dt.  (16) 

Jo 

The  exact  pdf,  peq ,  is  often  not  known  for  most  nonlinear  systems;  therefore, 
some  approximation  of  peq  is  needed  [23,  27,  25].  The  simplest  approximation 
is  to  use  the  Gaussian  pdf,  pfq(u,s),  with  the  same  mean  and  covariance  as 
in  the  original  system  [21,  14,  16,  2,  3,  4,  23,  27,  25].  This  approximation  is 
called  the  quasi-Gaussian  approximation  for  FDT  (qG-FDT).  The  correspond¬ 
ing  functional  Bp  becomes 


Bp{u,  s) 


V U{a(u,  s)pfq(u,  s)) 
P%(u,s) 


(17) 


and  the  the  quasi-Gaussian  approximation  to  the  response  function  (11)  is  given 

by 


R°{t)  =  ~(A(u(t  +  s),t  +  s)  <g>  Bp(u(s),s)j,  (18) 


while  the  corresponding  infinite  time  response  operator  becomes 

nOO 

1lG  =  RG(t)dt.  (19) 

Jo 

If  the  mean  and  the  variance  of  the  time-periodic  statistical  equilibrium  solution 
of  (3)  is  known  the  functional  in  (17)  can  be  computed  analytically  for  a  given 
perturbation  (9).  In  Section  4,  we  will  compute  the  functional  BG  for  the  test 
model  (1),  (2)  using  exactly  solvable  first  and  second  order  statistics  of  this 
model. 

Next,  we  present  a  few  examples  of  the  nonlinear  functionals  and  perturba¬ 
tions.  Very  general  nonlinear  functionals  have  the  separable  form 

A(u,s)  =  A(u)<j>(s),  (20) 
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where  <j>(s)  is  periodic  with  period  To.  For  the  response  of  the  mean  of  u,  we 
choose 

[A(u)]j=u3,  (21) 

while  for  the  response  of  the  variance,  we  take 

[j4(it)]j  =  ( Uj  —  ( Uj ))  .  (22) 

Substituting  (20)  into  (11)  and  using  the  change  of  variables  s'  =  s  +  t,  we 


rewrite  (11)  as 

1 

/•To 

II 

HO 

n. 

/  4>{t  +  s)(A(u(t  +  s))  <8>  Bf(u(s),  s)} 

Jo 

1 

pTo+t 

= 

n. 

/  <^(s,)(v4(u(s/))  <S>  Bp{u{s'  —  t),  s'  —  t))ds' 

by  periodicity 

1 

Tb. 

/•To 

/  </>(s')(A(u(s'))  <S>  Bp(u(s'  —  t),  s'  —  t))ds' .  (23) 

Jo 

Suppose  we  would  like  to  know  the  response  of  the  time  averaged  mean  or 
variance  over  a  part  of  the  period  [0,To],  say  over  [£i ,  ^2]  C  [0,To].  Practically, 
this  becomes  useful  when  monthly,  seasonal,  or  annual  averages  of  mean  or 
variance  are  of  interest.  Then,  the  choice  of  </>(s)  is  the  normalized  indicator 
function  of  the  segment  [ti ,  <^2] 

<t>{s)  =  ,  T°,  Xti,t3(g)>  (24) 

12  ~ 

where 

[  1,  for  s/mod(T0)  G  [ti,t2\, 

Xii,t2(s)=<  (25) 

I  0,  otherwise. 

Note  that  as  defined  Xti,t2  is  a  periodic  function  of  s  €  K1.  For  the  special 
choice  of  (f>(s)  given  in  (24),  the  response  function  in  (23)  becomes 
1  ft2 

R(t)  = -  /  {A(u(s))  ®  Bp{u(s —  t),s —  t))ds,  (26) 

h  ~  1 1  Jtl 

and  the  quasi-Gaussian  approximation  of  R(t)  is 

RG(t)  =  1  [  (A(u(s))  ®  Bp{u{s  —  t),s  —  t))ds.  (27) 

*  2  —  1 1  Jtl 
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Next,  we  discuss  two  types  of  perturbations  that  we  are  going  to  study  in 
this  paper:  perturbations  of  forcing  and  perturbations  of  linear  dissipation.  Per¬ 
turbations  of  forcing  are  described  by  (9)  with  a(u,  s)  =  1  and  the  components 
of  the  corresponding  functional  B f  becomes 

^UkPeqju,  S ) 

peq(u,s)  ■  [ 

On  the  other  hand,  the  perturbations  of  dissipation  are  given  by  a(u,  s)  =  —u 
with  the  functional  Bd  given  by  its  components 


[Bd]k 


^ uk  (^ kPeq  (n,  ^)) 
Peq(u,s) 


(29) 


In  the  quasi-Gaussian  approximation,  we  use  the  Gaussian  pdf,  pfq,  instead  of 
the  original  pdf,  peq 


and 


[Bf]k 


V UkPfq{u,S ) 
P%(U,S) 


[B% 


Vuk{Ukpfq{u,s)) 

P%(U,S) 


(30) 


(31) 


Now  following  [27],  we  discuss  practical  implementation  of  (11)  or  its  special 
case  (26)  for  a  given  dynamical  system  (3).  In  (26)  or  (27)  we  use  the  ergodicity 
of  system  (3)  and  substitute  the  phase  average  by  the  time  average  over  a  long 
time  trajectory  in  a  time-periodic  equilibrium  regime.  In  order  to  approximate 
the  probability  distribution  over  the  period  [0,Xo],  we  discretize  this  period 
with  L  equal  bins  centered  at  points  Sj.  Then,  the  response  function  (11)  can 
be  approximated  by 

1  l  „t*+t 

R(t)  w  —  S']  /  A{u{t  +  Sj  +  r),  t  +  Sj  +  t)  <gi  BF{u{sj  +  r),  Sj  +  t)cIt ,  (32) 
J-j-1  I T1* 
l=i  JT 

where  T  1,  and  T*  is  large  enough  so  that  the  system  (3)  has  reached 
time-dependent  equilibrium  at  this  time  T*.  Also,  we  always  use  L  =  12  cor¬ 
responding  to  a  monthly  partition.  Below,  we  will  test  the  quasi-Gaussian 
approximation  of  the  time-periodic  FDT  on  the  triad  model  (1),  (2). 
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It  is  worth  remarking  here  that  there  are  general  moment  constraints  in  the 
time-periodic  statistical  equilibrium.  These  constraints  allow  one  to  find  some 
higher  order  statistics  recursively  knowing  the  lower  order  statistics.  Then,  the 
response  of  these  higher  order  statistics  to  the  changes  in  forcing  or  dissipation 
can  be  found  automatically  once  the  response  of  the  lower  order  statistics  is 
found.  Consider  a  general  functional  G(u,  s)  which  is  periodic  in  s  with  the 
period  Tq.  Then,  we  find 


dG\ 

ds  / 

1 

: 


jr  J  J  -^(G(u,  s))peq(s)duds  = 
G(u,s)(s7(Fpeq)  -  -VV(Qpeq)^duds 


G(u,  s)—peq(s)dsdu 

~(fVg}  -  ^ct(7TVVG^ 
(33) 


To  illustrate  the  simplest  moment  constraints,  consider  G(u,s )  =  Go(u)  so 
that  (33)  becomes  the  moment  constraints 

(FVG0)  +  i(crerTVVGo)  =  0.  (34) 

Below  in  Eqs.  (56)-(59)  of  Section  3,  we  will  use  constraint  (34)  to  find  the 
annual-averaged  third  order  moment  that  has  an  important  physical  meaning 
for  the  system  (1),  (2)  —  it  represents  the  energy  transfer  among  the  modes. 


3.  Exactly  solvable  triad  model 

As  we  introduced  earlier  in  (1),  (2)  the  model  that  we  use  here  for  test¬ 
ing  time-periodic  FDT  is  a  stochastic  triad  nonlinear  model  with  time-periodic 
deterministic  forcing 

^  =  -71u1  +  fi{t)  +  CTiWi,  (35) 

^  =  (-72  +*(wo  +a0ui))u2  +  flit)  +  cr2W2,  (36) 

where  u\  is  real  mode  and  U2  is  a  complex  mode.  In  Eqs.  (35)  and  (36), 

the  parameters  71  and  72  represent  linear  dissipation  and  a  1  and  02  repre¬ 

sent  stochastic  forcing.  Also,  ojq  is  the  deterministic  linear  frequency  of  U2  and 


12 


ao  characterizes  the  strength  of  nonlinear  coupling  of  u 2  with  u\.  We  choose 
periodic  forcings  fi(t)  and  f2(t)  with  the  same  period  To  which  represents  sea¬ 
sonal  cycle.  As  we  discussed  in  the  introduction,  we  call  the  period  [0,  T0]  a 
“year”  as  in  climate  science.  Naturally,  we  divide  the  year  into  four  seasons 
with  each  season  consisting  of  three  months.  System  (35),  (36)  was  studied  by 
the  authors  in  [11,  12]  in  a  context  of  filtering.  There,  it  was  shown  that  there 
exist  analytical  formulas  for  the  first  and  second  order  statistics  of  tti  and  U2- 
Here,  we  will  use  those  formulas  to  compute  the  exact  “ideal”  response  of  the 
mean  and  the  variance  of  the  system  (35),  (36)  to  the  perturbations  of  forcing 
or  dissipation.  We  are  interested  in  the  infinite  time  response  of  the  system,  i.e. , 
we  measure  the  statistics  of  the  system  after  infinite  time  has  passed  since  the 
perturbation  was  applied.  Practically,  the  relaxation  time  which  is  much  longer 
than  the  system  decorrelation  time  has  to  pass  before  the  response  of  the  system 
is  measured.  Moreover,  the  time-periodic  equilibrium  mean  and  covariance  will 
be  used  to  construct  quasi- Gaussian  approximation  of  FDT.  The  ideal  response 
will  be  used  as  a  benchmark  for  the  qG-FDT  response. 

We  show  how  to  compute  the  mean  and  the  variance  of  tti  and  112  by  fol¬ 
lowing  [11,  12].  The  solution  of  Eq.  (35)  is  given  by 

ui(i)  =  uioe-7l(t-to)  +  [  /i(s)e-7l(t-s)ds  +  CT1  [  e"7l(t-s)dWi(s)(37) 

Jto  Jto 

where  itio  is  the  initial  condition  at  t  =  to-  Now  it  is  easy  to  find  the  mean  and 
the  variance  of  u\.  The  mean  is  given  by 

<«i(f)>  =  <«io)e-7l(‘-to)  +  f  fi(s)e~ll(-t~s^ds,  (38) 

Jto 

The  statistics  in  the  time-periodic  equilibrium  can  be  found  by  sending  to  to 
minus  infinity,  which  is  equivalent  to  allowing  an  infinitely  long  time  to  pass  after 
the  initial  condition  is  imposed.  Then,  the  time-periodic  equilibrium  mean  of 
Mi  is  given  by 

(ui  (t))eq  =  (  /i(s)e_7l(i_s)ds,  (39) 

J  —OO 

Similarly,  the  variance  of  u\  at  time  t  with  initial  condition  given  at  time  to 


13 


becomes 


Var(Ul(t))  =  y«r(wio)e-7l(*-to)  +  ^(1  -  e-7l(4-4o)). 

271 

The  equilibrium  variance  becomes 


.  .  CJ 1 

Vareq^Ux)  =  - — . 


(40) 


(41) 


eqK^'  27i' 

Note  that  ui(t)  is  Gaussian  and,  therefore,  it  is  fully  defined  by  its  mean  and 
variance. 

Now,  we  find  the  mean  and  the  covariance  of  the  second  mode,  u2-  The 
solution  of  Eq.  (36)  with  initial  condition  at  time  to  is  given  by 

ft  ft 


u2{t)=e  72(4  to)ip(t0,t)u2o  +  f  e  72(t  a)ip{s,t)f2{s)ds  +  (t2  f  e  72(4  ^ip(s,t)dW2{s),  (42) 

Jto  J  to 

where  as  in  [11,  12],  we  define  new  functions 

V>(s,f)  =  eiJ(s'4),  (43) 

J(s,t)  =  J  (coo  +  aoUi(s')^ds' =  (t,  —  s)uo  +  cio  J  Ui(s')ds'  (44) 

=  JD{s,t)  +  Jw(s,t) +  b(s,t)uw, 


where  the  deterministic  part  of  J(s,t)  is 

f4  rs' 

f  S  j  to 


JD{s,t)  =  (t-s)co0+a0  (  f  /i(s")e  7l(s  s  ]d s'ds ", 

J  S  Jto 

the  noisy  part  of  J(s,t)  is 

Jw(s,t)  =  criao  f  ds’  f  e7l(s""s'}(ilTi(s"), 

J  S  J  to 

and  the  prefactor  of  M10  is 

b(s,  t)  =  —( e-7l(s-4°)  -  . 

71  V  1 

The  mean  of  U2  is  given  by 

{u2(t))  =  e_72(4_to)(^(t0,t)w2o)  +  f  e~'l2<'t~s)(ip(s,t)}f2(s)ds. 

Jtn 


(45) 
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In  the  time-periodic  equilibrium  state  the  mean  of  u2  becomes 

{ui{t))eq  =  f  e~l2(t~s)(4>{s,t)}eqf2(s)ds.  (46) 

J  —  OO 

Note  that  J{s,t)  is  Gaussian  since  it  is  an  integral  of  the  Gaussian  variable  u\ 
plus  a  deterministic  function  as  can  be  seen  from  (44).  Therefore,  (ip(s,t))eq  is 
by  definition  a  characteristic  function  of  a  Gaussian  random  variable,  which  has 
the  explicit  form 

(4>{s,t))eq  =  {eiJ(s’t))eq  =  (47) 

The  mean  and  the  variance  of  J(s,  t)  were  computed  in  [11,  12]  and  here  we  use 
those  formulas 

(J(s,t))eq  =  (t  -  s)w0  +  a0  f  ds'  f  /i(s")e_7l(s,_s"W',  (48) 

J  s  J  —  OO 

and 

V areq(J (s,  t))  =  ^p(e_7l|t_s|  +  7i \t  ~  s\  -  l) ,  (49) 

where  again  we  used  to  — *  as  initial  time  to  find  time-periodic  equilibrium 

statistics.  Similarly,  we  can  find  Coveq(u2,ui),  Vareq(u2),  and  Coveq(u2,U2)- 
We  put  these  computations  in  the  Appendix. 

Although  expressions  (46)  and  (77),  (73),  (74)  are  analytical  and  explicit, 
the  integrals  still  have  to  be  computed  using  numerical  quadrature.  Exponential 
decay  of  the  integrands  in  (46),  (77),  (73),  (74)  suggests  that  we  change  variables 
in  order  to  replace  the  computation  of  the  integrals  over  a  semi-infinite  interval 
(— oo,t]  with  the  computation  of  the  integrals  over  a  finite  interval  [0, 1].  We 
define 


q  =  e-72(‘-s) 


which  gives 


s  =  t  H - log(g). 

72 


(50) 


(51) 
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(52) 


Then,  Eq.  (46)  becomes 

1  r1 

(u2(t))eq  =  —  /  ( 4>{s(q),t))eqf2(s(q))dq . 

72  Jo 

The  same  approach  is  used  in  computing  the  cross-covariance  Cov(u2,  u\)  given 
by  (77).  Similarly,  we  deal  with  the  double  integrals  for  the  second  order  statis¬ 
tics  in  (73),  and  (74).  The  change  of  variables  (51)  is  equivalent  to  using  a 
non-uniform  (logarithmic)  mesh  in  the  quadrature  for  the  original  variables. 
Next,  we  apply  the  trapezoidal  rule.  We  partition  the  interval  [0, 1]  with  points 
qj  such  that 

0  =  qo  <  qi  <  ■  ■  ■  <  qx  =  1,  (53) 


with  a  small  uniform  step  A q  =  qj+i  —  qj.  Then,  the  integral  of  any  function 
g(q)  over  the  interval  [0, 1]  is  approximated  according  to 

,1  ^  K—l 

/  g(q)dq  «  -y-(ff(<?o)  +9(qk))  +  A<?  9(Qj)-  (54) 

J0  i=i 

Note  that  in  all  our  examples,  the  integrand  g{q)  vanishes  for  q  =  0  because 
of  the  exponential  decay  to  zero  of  (ip(s,t)}  as  s  — >  — oo  (which  is  equivalent 
to  q  — *  0).  It  is  convenient  to  introduce  real  variables  of  the  triad  system  (35) 
and  (36) 


( X1) 

f  U1  \ 

u  = 

X2 

= 

Re[u2] 

V  *3 ) 

\  Im[u2]  y 

(55) 


In  Fig.  1,  we  compare  the  first  and  second  order  statistics  of  the  triad  (55), 
computed  in  this  Section  analytically  with  the  corresponding  results  of  Monte- 
Carlo  ensemble  averaging.  The  conversion  formulas  for  the  statistics  of  the  real 
triad  (55)  from  the  statistics  of  the  complex  system  {u\,u2)  are  given  in  the 
Appendix.  Here,  we  used  the  parameter  set  given  in  the  second  row  of  Table  1 
(this  Table  will  be  discussed  below  in  Section  5).  The  Monte-Carlo  averaging 
was  done  using  a  very  long  trajectory  (36  years)  computed  via  Eqs.  (37),  (42). 
We  used  a  time  step  h  =  10-3  to  approximate  the  integrals  in  (46)  and  in  (77), 
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(73),  (74).  Note  that  throughout  paper  the  time  is  measured  in  either  years  or 
months  here,  where  1  year  is  equal  to  2.4  time  units  in  terms  of  the  coefficients 
given  in  Table  1  and  1  month  is  1/1 2-th  of  the  year. 

We  note  excellent  agreement  between  the  analytical  formulas  and  the  results 
of  Monte-Carlo  averaging  shown  in  Fig.  1.  Moreover,  in  Fig.  1,  we  demonstrate 
the  skewness  and  flatness  of  the  nonlinear  variable  X2  computed  using  Monte- 
Carlo  averaging.  We  note  a  large  burst  of  non-Gaussianity  in  the  third  season 
of  the  annual  cycle. 

Next,  we  apply  Eq.  (34)  to  compute  the  annual-averaged  third  order  mo¬ 
ment  (272:22:3)  that  represents  the  simultaneous  energy  transfer  among  the  three 
modes  of  the  system  (35),  (36)  by  utilizing  combination  of  the  first  and  second 
order  moments.  For  any  functional  Gq(u),  Eq.  (34)  becomes 

+  ((“72^2  ^  W0X3  -  CI0X1X3  +  Re[f2\{t))^^j 

+({-12X3  +u>0x2  +a0xiX2  + 

+^(diag[a?,<ri/2,®2/2]Vu  •  VuG0(u,i))  =  0,  (56) 

where  diag[-]  is  a  diagonal  matrix  with  the  given  elements.  Let  us  apply  this 
equation  to  some  specific  functionals.  For  G 0  =  £§  we  find 

(27272:3)  =  —  (72(2:3)  -  u}0\x2X3)  -  (/m[/2](t)x3)  -  ^2.) .  (57) 

CLq  \  4  / 

Similarly,  for  G  =  x\  we  have 

(2727x3)  =  —  (  -  72(2:2)  -  w0(x2x3)  +  \Re[f2}{t)x2)  +  ^2.)  ■  (58) 

CLq  \  4  / 

We  add  these  two  expressions  to  find 

(272:22:3)  =  ^(  -  y^)  -  (2:2))  -  w0(x2x3)  +  i((i?e[/2](t)x2)  -  (/TO[/2](t)x3)).  (59) 

Note  that,  the  right  hand  side  of  Eq.  (59)  contains  only  the  first  and  the  second 
order  statistics  which  can  be  computed  as  discussed  above.  The  annual-averaged 
triple-correlator  (x  1X2X3)  is  an  important  physical  quantity  that  measures  the 
energy  exchange  among  all  three  modes  simultaneously,  which  characterizes  the 
nonlinear  wave  coupling. 
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The  analytical  formulas  for  the  first  and  second  order  statistics  can  be  ap¬ 
plied  to  the  computation  of  the  ideal  response  of  the  triad  system  to  the  external 
perturbation.  In  particular,  in  this  paper  we  will  be  studying  the  mean  response 
to  the  changes  in  forcing  and  the  variance  response  to  the  changes  in  dissipation. 
The  ideal  response  is  the  actual  linear  response  to  the  perturbation  and  here 
we  know  it  explicitly  for  the  mean  and  variance  unlike  the  extensive  numerical 
calcualtion  for  this  response  required  for  other  systems  ([23,  22,  14,  16]).  The 
ideal  response  operator  computed  below  will  be  used  as  a  benchmark  for  the 
corresponding  qG-FDT  response  operators  that  are  discussed  in  Section  4.  We 
start  with  the  ideal  mean  response  to  the  change  in  forcing.  For  simplicity,  we 
drop  the  subscript  “eq”  in  our  notation  assuming  that  all  averages  are  done  in 
the  time-periodic  statistical  equilibrium.  We  define  the  real  forcing  vector 


/  = 


(  f  \ 

J  X  i 

fx  2 


\  fx  3  / 


(60) 


The  ideal  mean  response  to  the  changes  in  forcing  is  give  by  the  matrix 


d(xi) 

d(xi) 

d(xi) 

9fxK 

d/x2_ 

d/x3_ 

d~(x2 ) 

d~(x2) 

d~(x2) 

9fxi_ 

d_fx  2_ 

d/*3_ 

d~(x3j 

d~{x3j 

d~(x3) 

9fxi 

dfx2 

dfx3 

(61) 
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Using  Eqs.  (39)  and  (46),  we  find 


1 

7 

7i 

=  0, 

=  0, 

= - ,  a°  Jm  [  [  e~'l2{t~s) (if{s,t)){t  -  s)f2{s)dsdt, 

71U2-U)  Jtl  7-oo 

=  ,  a°  ,  Re  [  [  e~l2it~s) -  s)f2{s)dsdt, 

71U2-U)  Jtl  7-00 

=  ^rL  =  T^TRef  2  f  ^-^WMdsdt, 

OJx 3  t2  ~  t\  Jtl  J_ob 

=  -dJf3  =  - .  Im  [  [  e_72(t_s)(^(s,t))dsdi. 

O/aa  12  -  El  Jtl  7-oo 

In  the  Appendix,  we  show  how  to  find  the  ideal  operator  for  the  variance  re¬ 
sponse  to  the  changes  in  external  forcing,  TZyj,  and  dissipation,  7 Zyd. 

4.  Quasi-Gaussian  approximation  of  FDT  for  the  triad  system 

In  Section  2,  we  briefly  provided  a  general  theory  for  the  time-periodic  FDT. 
There,  we  obtained  the  expression  for  the  FDT  response  operator  given  by 
Eq.  (11).  However,  if  the  time-periodic  equilibrium  pdf,  peq,  is  not  known 
as  happens  in  most  practical  approximations,  the  quasi-Gaussian  approxima¬ 
tion  (18)  to  the  response  operator  is  computed  using  the  Gaussian  pdf,  pfq  with 
the  same  mean  and  covariance  as  the  original  system. 

Here,  we  compute  the  qG-FDT  response  operator,  1ZG ,  for  the  triad  model  (55) 
for  the  time-averaged  mean  and  variance  response  to  the  changes  of  forcing  or 
dissipation.  We  will  be  using  different  averaging  windows  that  correspond  to 
the  first  month,  the  third  season,  and  the  full  year  with  formulas  as  developed 
in  (24)- (27). 

For  the  mean  response,  the  corresponding  functional  A(u)  in  (27)  has  the 
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dfXl 
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form 


Am{u) 


and  for  the  variance  response 


(Xl\ 

X2 

\x3  ) 


Av{u) 


^  (xt  -  (xi))2  ^ 
(X2  -  (X2))2 

y  (x3  ~  (^3))2  ) 


(62) 


(63) 


Now,  we  compute  the  functional  Bp  using  (30)  and  (31)  for  the  perturbations 
of  forcing  and  dissipation,  respectively.  Note  that  the  Gaussian  pdf,  pfq,  with 
the  same  mean  and  variance  as  the  solution  of  the  original  system  (35)  and  (36) 
in  time-periodic  equilibrium  is  computed  analytically  using  analytical  formulas 
for  the  first  and  second  order  statistics  of  (35)  and  (36)  obtained  in  Section  3 


PG(u,s)  =  ^3|S|  exp  ^-^(u-  (m))tS  ^-(w})^,  (64) 

where,  the  mean,  (it),  and  the  covariance  matrix,  S,  of  the  triad  system  are 
computed  in  the  Appendix.  For  the  perturbations  of  forcing  we  use  (30)  and  (64) 
to  find 


[Bf(u)]k  =  yk  =  \B  1(u-  (u))]fc. 


(65) 


Then,  the  quasi-Gaussian  approximation  (18)  of  the  mean  response  function  to 
the  changes  of  forcing  is  given  by  the  matrix  with  the  elements 


rt  2 


(■ AM)j{u(s))Bfk(u(s  —  t),s  —  t))ds 


h  - 1\  Jtl 

i  rt2 

— T  /  (xj(s)yk(s  —  t))ds. 

*1  JU 


h  -  <u  Jt ! 

Next,  we  consider  perturbations  of  dissipation  by  combining  (31)  and  (64) 
[B$(u)\k  =  1  -  £fc[£-1(w  -  (it))]fc  =  1  -  XkVk- 


(66) 


(67) 
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Then,  the  quasi-Gaussian  approximation  (18)  of  the  mean  response  operator  to 
the  changes  of  dissipation  is  given  by  the  matrix  with  the  elements 

[RM,d\jk(t)  =  1  f  (xj(s)(l  -  Xk{s  -  t)yk{s  -  t)))ds.  (68) 

1 2  -  1 1  Jtl 

Similarly,  we  find  the  variance  response  to  changes  in  forcing 

iRvj]jk(t)  =  1  [  (( Xj(s )  -  (xj(s)))2yk(s  -  t))ds,  (69) 

1 2  -  tl  Jtl 

and  in  dissipation 

[Rv,d\jk(t)  =  ,  1  ,  /  ((®j(s)-(a:j(s)))2(l-a;fc(s-t)j/fc(s-t)))ds.  (70) 

C2  —  ti  ytl 

The  infinite  time  response  operator,  77° ,  is  computed  from  the  response  func¬ 
tions,  i?G(t),  via  Eq.  (19).  We  note  that  the  quasi-Gaussian  response  function 
for  the  mean  to  the  change  of  forcing,  Rfj  ^(t),  is  basically  the  second  order  two- 
time  correlation  function,  for  the  mean  to  the  changes  in  dissipation,  Ryj(t), 
and  the  variance  to  the  changes  in  forcing,  d(t),  are  the  third  order  two-time 
correlation  function,  and  the  variance  to  the  changes  in  dissipation,  Ryd{t),  is 
the  fourth  order  two-point  correlation  function.  Moreover,  the  odd-order  cen¬ 
tered  statistics  of  Gaussian  random  variables  vanish  [10],  and  hence,  for  the 
near-Gaussian  regimes,  the  mean  response  to  the  perturbations  of  dissipation 
and  variance  response  to  the  perturbations  of  forcing  do  not  deviate  from  zero 
significantly.  Therefore,  in  our  study  we  will  concentrate  on  the  mean  response 
to  the  perturbations  of  forcing  and  variance  response  to  the  perturbations  of 
dissipation. 

If  we  assume  that  the  perturbations  of  forcing,  Sf,  or  perturbations  of  dissi¬ 
pation,  Sj,  have  three  independent  components,  Sfj  or  Sjj,  that  correspond  to 
each  of  the  variables  Xj  then  the  corresponding  response  operators,  7 Z  are  3x3 
constant  coefficient  matrices.  However,  while  this  assumption  of  independence 
of  perturbation  components  is  valid  for  the  perturbations  of  forcing,  it  is  not 
true  for  the  perturbations  of  dissipation  because  we  only  have  two  dissipation 
parameters,  71  and  72.  As  we  show  in  the  Appendix,  the  variance  response  to 
the  changes  in  dissipation  is  described  by  a  3  x  2  constant  coefficient  matrix. 
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In  order  to  find  the  corresponding  qG-FDT  response  operator,  we  equate  the 
perturbations  of  dissipation  for  the  second  and  the  third  components  of  the 
triad  (55),  X2  and  X3.  This  is  equivalent  to  adding  the  second  and  the  third 
columns  of  the  qG-FDT  response  matrix  that  is  obtained  under  the  assump¬ 
tion  of  independence  of  all  three  perturbation  components.  Then,  the  qG-FDT 
variance  response  to  the  changes  in  dissipation  is  described  by  a  3  x  2  constant 
coefficient  matrix. 

5.  Skill  of  the  quasi-Gaussian  FDT  on  the  test  model 

In  this  Section,  we  study  the  skill  of  the  qG-FDT  response  for  the  triad 
system  (35),  (36)  in  comparison  with  the  ideal  response.  Here,  the  terminology 
skill  refers  to  the  capability  of  the  quasi-Gaussian  FDT  approximation  to  re¬ 
produce  the  exact  ideal  response  reported  for  example  below  (61)  for  the  mean. 
We  have  chosen  a  number  of  different  parameter  regimes  presented  in  Table  1 
for  our  tests.  As  the  first  test-case,  we  consider  a  near-Gaussian  regime  which 
is  the  most  obvious  for  testing  qG-FDT  because  the  qG-FDT  provides  an  exact 
ideal  response  when  it  is  applied  to  Gaussian  systems.  As  the  second  test-case, 
we  used  the  regime  when  the  system  (35)  and  (36)  has  near-Gaussian  statistics 
over  the  whole  year  except  for  the  third  season,  when  the  pdf  deviates  from 
Gaussian  significantly.  In  this  regime  we  expect  deteriorating  skill  for  qG-FDT 
when  the  time  averaging  of  the  statistics  is  done  over  this  non-Gaussian  season. 
Otherwise,  we  still  anticipate  high  skill  of  qG-FDT.  In  the  first  two  cases,  we 
compared  the  qG-FDT  and  ideal  mean  responses  to  the  changes  of  forcing.  For 
the  third  test-case,  we  use  a  series  of  simulations  for  the  triad  system  with  a 
systematically  increasing  nonlinearity  with  a  forcing  that  is  non-resonant.  In 
this  case,  we  expect  qG-FDT  to  lose  skill  as  the  nonlinearity  increases  and  the 
system  becomes  more  non-Gaussian.  Finally,  the  fourth  test  case  is  the  same 
as  the  third  one  but  with  the  resonant  forcing.  Here,  the  resonant  forcing  am¬ 
plifies  the  statistics  of  the  model  and  hence  the  model  has  large  variance  and 
becomes  strongly  non-Gaussian  and,  therefore,  the  qG-FDT  is  not  expected  to 
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have  much  skill  as  nonlinearity  increases.  In  the  third  and  fourth  test  cases  we 
compare  the  mean  response  to  the  changes  of  forcing  and  the  variance  response 
to  the  changes  of  dissipation.  Below  we  provide  a  detailed  explanation  of  the 
testing  procedure  presented  for  the  first  test  case  with  the  near-Gaussian  regime. 
For  the  remaining  test  cases  we  use  the  same  methodology  and  only  discuss  the 
results  of  our  study. 

5.1.  qG-FDT  response  in  the  near-Gaussian  regime 

The  first  test-case  represented  in  Table  1  is  characterized  by  weak  damping, 
weak  white  noise  forcing  and  weak  deterministic  forcing  of  the  nonlinear  mode, 
U2-  This  regime  is  very  close  to  Gaussian  and  the  quasi-Gaussian  approximation 
to  FDT  is  expected  to  perform  well  here. 

We  compare  the  mean  qG-FDT  response  to  the  changes  in  forcing  with 
the  corresponding  ideal  response  in  this  near-Gaussian  regime.  As  discussed  in 
Section  3,  we  utilize  the  analytical  formulas  (39)  and  (46)  of  the  time-periodic 
equilibrium  mean  to  find  the  ideal  mean  response  of  the  system  (35),  (36)  to  the 
perturbations  of  forcing.  We  use  three  types  of  time  averaging,  i.e.,  monthly 
average  over  the  first  month  (ti  =  0  and  1 2  =  To/12  in  Eq.  (24)),  seasonal 
average  over  the  third  season  (t\  =  T0/2  and  1 2  =  3T0/4  in  Eq.  (24)),  and 
annual  average  (t\  =  0  and  £2  =  T0  in  Eq.  (24)).  We  construct  the  response 
matrix  following  the  procedure  outlined  in  Section  3.  These  ideal  response 
operators  along  with  their  singular  values  are  given  in  the  upper  left  quarter  of 
Table  2. 

Next,  we  test  how  well  the  qG-FDT  predicts  the  ideal  linear  response. 
The  quasi-Gaussian  approximation  of  the  response  operator,  1ZG ,  is  given  by 
Eq.  (19).  First,  we  have  to  compute  the  response  function,  Rclfj{i),  which 
is  given  by  Eq.  (18)  with  the  functional  A  from  (62)  and  the  functional  B 
from  (30).  The  response  function,  f{t),  is  a  two-time  correlation  function 
and  we  use  the  ergodicity  of  system  (35),  (36)  to  compute  RGj  f(t)  averaging 
over  a  long  trajectory.  We  apply  the  computational  algorithm  “on  the  fly”  as 
described  in  [23]  when  the  correlation  function  is  computed  in  parallel  with  to 
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the  computation  of  the  trajectory  and  thus  a  very  long  trajectory  can  be  used 
in  order  to  achieve  a  very  high  precision. 

In  Fig.  2,  we  show,  Rfrj(t).  for  each  of  the  components  of  the  operator  for 
the  mean  to  the  changes  in  forcing.  Here,  we  used  a  long  trajectory  of  the  length 
3  •  106  oscillation  periods  (years).  Next,  following  (19)  we  integrate,  j(t) 
and  obtain  the  components  of  the  linear  response  matrix,  TZ^  ^  along  with  the 
corresponding  singular  values  given  in  the  upper  right  quarter  of  Table  2. 

We  compare  the  ideal  response  matrix  with  the  qG-FDT  response  matrix  to 
find  excellent  agreement  as  we  expected  for  this  near-Gaussian  regime.  More¬ 
over,  we  note  that  annually  averaged  operators  have  the  best  agreement  while 
monthly  averaged  operators  have  the  worst  agreement  among  the  three  types 
of  time  averaging.  This  can  be  explained  by  the  fact  that  averaging  over  one 
certain  month  requires  a  much  longer  trajectory  than  averaging  over  the  whole 
year  for  the  same  precision.  To  summarize,  we  have  confirmed  that  the  qG-FDT 
gives  excellent  prediction  of  the  mean  response  operator  for  the  perturbations 
of  forcing  for  the  system  in  a  near-Gaussian  regime. 

5.2.  qG-FDT  response  in  the  near-Gaussian  regime  with  non-Gaussian  season 

The  second  test-case  represented  in  Table  1  is  characterized  by  stronger 
damping  of  both  modes  and  stronger  white  noise  forcing  of  the  first  mode  when 
compared  with  the  first  test-case  from  Table  1.  Moreover,  in  the  second  test- 
case  we  have  stronger  forcing  / 2  of  the  nonlinear  mode  U2  which  creates  stronger 
non-gaussianity  in  the  statistics  of  U2 ■  In  Fig.  1,  which  we  already  studied  in 
Section  3,  we  demonstrate  the  statistics  of  the  triad  (55)  for  this  test-case. 
In  Section  3,  we  used  Fig.  1  to  confirm  that  the  analytically  obtained  first 
and  second  order  statistics  coincide  with  the  corresponding  results  of  Monte- 
Carlo  averaging.  Moreover,  the  skewness  plot  in  Fig.  1  shows  that  the  system 
is  in  the  near-Gaussian  regime  for  most  of  the  time  except  for  the  burst  of 
non-Gaussianity  in  the  third  season  (between  months  6  and  9).  In  Fig.  3,  we 
show  the  pdfs  of  X2  for  our  second  test-case  averaged  over  the  first  month, 
the  third  season  and  the  whole  year.  These  pdfs  were  computed  in  the  time- 
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periodic  equilibrium  and  rescaled  to  have  mean  zero  and  variance  one.  Note 
that  the  pdfs  here  and  throughout  the  rest  of  the  paper  are  computed  using 
Monte  Carlo  binning  with  80  bins  at  each  time.  The  pdfs  here  are  used  for 
qualitative  demonstration  the  non-Gaussian  properties  of  the  system,  on  the 
other  hand,  the  quantitative  analysis  of  the  behavior  of  the  qG-FDT  algorithm  is 
performed  using  the  exact  analytical  formulas  for  the  ideal  response  as  discussed 
in  Section  3.  We  note  that  the  pdf  of  X2  averaged  over  the  first  month  and  over 
the  whole  year  are  very  close  to  Gaussian  with  the  corresponding  skewness  —0.07 
and  0.05,  respectively.  On  the  other  hand,  the  pdf  of  X2  averaged  over  the  third 
season  is  significantly  non-Gaussian  with  the  skewness  0.57.  The  source  of  the 
non-Gaussianity  in  the  pdf  for  X2  can  be  traced  to  the  triplr  correlation  for 
anomalies,  x'j  =  Xj  —  (xj)  for  j  =  1,2,3,  i.e. ,  {x^x^x^};  this  triple  correlation 
vanishes  identically  for  a  Guassian  random  field.  In  Fig.  1,  we  plot  the  skewness 
of  X2  and  the  triple  correlator,  (x^x^x^),  as  functions  of  time  over  one  period. 
As  the  reader  can  see  by  inspection,  these  two  functionals  are  highly  correlated 
with  pattern  correlation  0.91.  In  the  Appendix,  we  present  the  analytic  formula 
for  (xix^x^). 

Now,  we  compare  the  mean  response  to  the  perturbations  of  forcing,  com¬ 
puted  using  ideal  and  qG-FDT  algorithms  as  discussed  in  Section  5.1.  We  refer 
to  the  lower  part  of  Table  2  now.  On  the  left  side,  we  see  the  ideal  response 
operator  together  with  the  corresponding  singular  values.  Again,  we  used  the 
averaging  over  one  month  (the  region,  bounded  by  the  dashed  vertical  lines  in 
the  second  panel  in  Fig.  1),  over  one  season  (the  region,  bounded  by  the  dotted 
vertical  lines  in  the  second  panel  in  Fig.  1),  and  over  one  year  (the  whole  pe¬ 
riod).  On  the  right  side  of  the  lower  part  of  the  Table  2,  we  show  the  qG-FDT 
mean  response  operators  to  the  changes  of  forcing  along  with  the  correspond¬ 
ing  singular  values.  We  note  that  qG-FDT  response  operator  approximates 
the  ideal  response  operator  quite  well  when  monthly  and  annual  averages  are 
considered.  As  we  mentioned  above  along  these  averaging  segments,  the  triad 
system’s  statistics  are  near-Gaussian.  However,  the  qG-FDT  operator  averaged 
over  the  third  season  is  significantly  different  from  the  corresponding  ideal  op- 
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erator.  In  particular,  we  note  the  difference  between  the  two  smallest  singular 
values  shown  in  bold  in  Table  2.  This  discrepancy  is  a  consequence  of  strong 
non-Gaussianity  of  the  model  in  the  third  season. 

To  conclude  the  results  for  the  second  test-case,  we  have  found  that  for  the 
system  with  near-Gaussian  statistics,  the  qG-FDT  yields  a  very  good  approxi¬ 
mation  of  the  ideal  response  operator,  while  in  the  non-Gaussian  regime,  there 
is  a  discrepancy  between  the  ideal  and  qG-FDT  mean  response  to  the  changes 
in  forcing. 

5.3.  qG-FDT  response  with  the  near-resonant  forcing 

In  this  and  the  next  Sections  we  study  how  the  skill  of  the  quasi-Gaussian 
approximation  of  FDT  changes  as  the  system  undergoes  a  transition  from  a 
Gaussian  to  a  non-Gaussian  regime  as  the  nonlinearity  becomes  stronger.  In 
this  Section,  we  consider  non-resonant  regimes  when  the  frequency  of  the  forcing 
uj  =  2tt/Tq  is  different  from  the  frequency  wo  of  the  nonlinear  mode  112-  We 
gradually  increase  the  strength  of  the  nonlinear  coupling  which  is  controlled  by 
the  parameter  a 0  and  study  how  much  the  skill  of  qG-FDT  changes.  Here,  we 
consider  not  only  the  mean  response  to  the  changes  in  forcing  given  by  (19),  (66) 
but  also  the  variance  response  to  the  changes  in  dissipation  given  by  (19),  (70). 
The  parameters  for  the  third  test-case  are  given  in  the  third  row  of  Table  1. 

In  Fig.  4,  we  show  the  evolution  of  the  pdfs  for  the  non-resonant  case  for 
a 0  =  0.5.  The  pdfs  are  shown  at  the  beginning  of  each  of  the  twelve  months.  We 
note  that  the  pdfs  here  are  near-Gaussian  by  sight.  We  also  note  the  floating 
mean  position  roughly  coincides  with  the  position  of  the  peaks  of  the  pdfs.  On 
the  other  hand,  the  variance  does  not  change  much.  In  order  to  quantify  non- 
gaussianity,  in  Table  3,  we  show  the  dependence  of  skewness  and  flatness  of  the 
triad  model  in  a  time-periodic  equilibrium  state  as  the  nonlinearity  parameter 
a 0  increases.  We  used  monthly,  seasonally,  and  annually  averaged  pdfs  for 
computing  skewness  and  flatness.  We  note  that  in  the  regime  with  non-resonant 
forcing,  the  system  stays  very  close  to  Gaussian  even  for  as  high  values  of 
nonlinear  parameter  as  ao  =  0.5.  Therefore,  the  qG-FDT  should  have  high 
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skill  here.  However,  as  we  will  see  in  Section  5.4  the  situation  changes  with 
resonant  forcing  when  the  system  becomes  strongly  non-Gaussian  even  with 
small  values  of  the  nonlinearity  parameter  a o- 

To  test  the  skill  of  the  qG-FDT,  we  followed  the  procedure  outlined  in  Sec¬ 
tion  5.1,  where  we  compared  the  ideal  and  qG-FDT  response  operators  and 
their  singular  values.  Here,  we  only  compare  the  singular  values  of  the  corre¬ 
sponding  operators.  In  Table  4,  we  demonstrate  the  singular  values  of  the  ideal 
and  qG-FDT  operators  for  the  time-averaged  mean  response  to  the  changes 
in  forcing.  We  note  that  the  corresponding  singular  values  are  very  close  to 
each  other  for  all  values  of  a o-  Moreover,  the  qG-FDT  approximations  to  the 
response  operators  for  the  monthly  averaged  mean  have  less  skill  than  the  re¬ 
sponse  operators  averaged  over  one  season,  which  in  turn  have  less  skill  than 
the  operators  for  the  annually  averaged  mean.  As  we  mentioned  in  Section  5.1, 
this  is  a  consequence  of  the  fact  that  we  use  the  same  trajectory  to  compute  the 
response  operators  for  all  three  time  averaged  statistics,  however,  the  averaging 
over  one  month  requires  more  data  to  achieve  the  same  precision  in  the  corre¬ 
lation  function  than  the  average  over  one  season  or  one  year.  Next,  we  consider 
the  variance  response  to  the  changes  in  dissipation.  We  compute  the  ideal  vari¬ 
ance  response  using  Eqs.  (83)  from  the  Appendix.  As  discussed  in  Section  3, 
the  operator  for  the  variance  response  to  the  changes  in  dissipation  is  of  the 
size  3x2  and,  therefore,  it  has  2  singular  values.  When  studying  the  variance 
response  to  the  changes  in  dissipation,  we  should  keep  in  mind  that  the  variance 
in  the  non-resonant  case  is  almost  independent  of  time  with  just  a  small  devia¬ 
tions  around  its  annually  averaged  value.  These  annual-averaged  values  of  the 
variance  are  given  in  Table  3  for  different  values  of  nonlinearity.  We  note  that 
the  annual-averaged  variance  grows  only  a  little  as  ao  increases.  Below,  in  Ta¬ 
ble  6  for  the  resonant  forcing  regime,  we  will  observe  a  much  more  rapid  growth 
of  the  annual-averaged  variance  as  a  function  of  a q.  In  Table  5,  we  compare 
the  singular  values  of  the  ideal  and  qG-FDT  variance  response  operators  to  the 
changes  in  dissipation  for  different  values  of  nonlinearity  parameter  ao  for  the 
case  of  non-resonant  forcing.  From  Table  5,  we  learn  that  the  skill  qG-FDT  ap- 
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proximation  of  the  variance  response  to  the  changes  in  dissipation  deteriorates 
as  the  strength  of  nonlinearity  increases.  In  particular,  we  note  that  the  larger 
singular  value  is  predicted  with  reasonable  accuracy  (with  ~  5%  error)  only  for 
the  nonlinear  parameter  values  up  to  do  =  0.3.  Then,  as  do  increases  further, 
the  larger  singular  value  of  the  qG-FDT  operator  is  significantly  larger  than  that 
of  the  ideal  operator.  This  discrepancy  for  large  nonlinearity  can  be  attributed 
to  the  fact  that  the  variance  response  function  (70)  is  the  fourth  order  two-time 
correlation  function,  whereas,  the  quasi-Gaussian  approximation  is  designed  to 
fit  the  first  and  the  second  one-time  moments.  On  the  other  hand,  the  smaller 
singular  value  is  predicted  quite  well  for  all  values  of  do- 

To  summarize,  we  have  tested  the  qG-FDT  approximation  of  the  mean  re¬ 
sponse  to  the  changes  in  forcing  and  the  variance  response  to  the  changes  in 
dissipation  on  the  non-resonant  forcing  test-case.  We  have  seen  that  the  mean 
response  to  the  changes  in  forcing  is  predicted  very  well  for  a  wide  range  of 
nonlinearity  parameter  do,  whereas  the  variance  response  to  the  changes  in  dis¬ 
sipation  slightly  deteriorates  when  the  nonlinearity  parameter  exceeds  certain 
value. 

5-4-  qG-FDT  response  with  the  resonant  forcing 

Now,  we  study  the  test-case  of  resonant  forcing  with  the  parameters  given 
in  the  fourth  row  of  Table  1.  In  Fig.  5,  we  demonstrate  the  monthly  snapshots 
of  the  evolution  of  the  time-periodic  pdf,  p(x 2 ,  s) ,  of  the  nonlinear  variable  X2 
of  the  triad  u  given  by  (55).  Figure  5  was  produced  in  the  strongly  nonlinear 
regime  with  do  =  0.35.  We  note  that  the  pdf  in  Fig.  5  has  two  peaks  during 
parts  of  the  annual  cycle.  In  Fig.  6,  we  show  the  time-averaged  pdfs  rescaled 
to  have  mean  zero  and  variance  one  for  increasing  values  of  do  and  for  three 
averaging  periods  over  the  first  month,  the  third  season,  and  the  full  year.  We 
compare  these  pdfs  with  the  standard  Gaussian  distribution  N{ 0, 1).  We  note 
that  monthly  averages  are  more  non-Gaussian  than  seasonal  averages  which,  in 
turn,  are  more  non-Gaussian  than  annual  averages.  To  monitor  the  departure 
from  Gaussian  statistics,  we  measure  corresponding  skewness  and  flatness  of 
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X2  and  present  them  in  Table  6.  Strong  deviations  of  the  system’s  statistics 
from  Gaussian  values  show  that  this  is  a  very  tough  test  case  for  the  qG-FDT. 
Moreover,  in  Table  6,  we  show  the  annual  averages  of  the  variance  of  x 2  as 
nonlinearity  increases.  We  note  a  much  faster  growth  of  the  annually-averaged 
variance  as  a  function  of  ao  than  it  was  in  the  case  with  non-resonance  forcing 
(Table  3).  The  very  rapid  growth  of  the  variance  of  X2  as  a  function  of  the 
nonlinearity  strength  ao  is  a  consequence  of  resonant  forcing. 

Next,  we  study  the  skill  of  the  qG-FDT  response  operators.  Similar  to  the 
non-resonant  case  discussed  in  Section  5.3,  we  first  discuss  the  mean  response 
to  the  changes  in  forcing.  In  Table  7,  we  show  the  singular  values  of  the  ideal 
and  qG-FDT  response  operators  for  three  averaging  windows  (first  month,  third 
season,  and  full  year)  and  for  increasing  strength  of  the  nonlinearity  parame¬ 
ter  a o-  We  note  that  as  the  nonlinearity  grows,  the  largest  singular  value  grows 
significantly  for  the  monthly  and  seasonal  ideal  averaged  mean  responses.  More¬ 
over,  the  other  two  singular  values  also  depend  on  the  nonlinearity  parameter 
a 0  for  the  monthly  and  seasonal  ideal  averaged  mean  response  operators,  one  of 
them  decreases  and  the  other  one  increases  as  ao  grows.  On  the  other  hand,  the 
singular  values  of  the  operator  for  the  annually  averaged  mean  response  to  the 
perturbations  of  forcing  are  practically  independent  of  ao .  This  phenomenon 
can  be  attributed  to  the  fact  that  annually-averaged  pdfs  of  the  system  are 
much  more  Gaussian  than  the  monthly  or  seasonally  averaged  pdfs  as  can  be 
seen  from  Fig.  6  and  Table  6.  In  Table  7,  we  compare  the  ideal  and  qG-FDT 
operators  for  the  time-averaged  mean  response  to  the  changes  in  forcing  and  the 
corresponding  singular  values.  We  also  show  these  singular  values  as  function  of 
nonlinearity  strength  ao  in  Fig.  7  (upper  panel).  We  note  surprisingly  excellent 
agreement  between  the  ideal  and  qG-FDT  response  operators  when  compared 
by  the  corresponding  singular  values.  Even  with  the  strongest  nonlinearity  that 
we  considered  (a 0  =  0.35),  the  error  in  predicting  the  singular  values  by  qG- 
FDT  does  not  exceed  10%.  This  numerical  experiment  shows  very  high  skill  of 
the  quasi-Gaussian  approximation  of  FDT  for  the  mean  response  to  the  changes 
in  forcing  even  in  a  strongly  non-Gaussian  regime  with  the  two-peaked  pdf. 
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Next,  we  consider  the  variance  response  to  the  changes  in  dissipation.  We 
found  that  the  variance  stays  almost  constant  over  the  whole  year  period,  there¬ 
fore,  the  annually  averaged  values  of  variance  given  in  Table  6  represent  the 
variance  in  the  each  month  of  the  year.  We  note  that  the  variance  increases  six 
times  as  a o  grows  from  ao  =  0  to  ao  =  0.35,  which  is  a  consequence  of  the  reso¬ 
nant  forcing.  In  Table  8,  we  compare  the  ideal  and  qG-FDT  variance  response 
operators  to  the  changes  in  dissipation  and  corresponding  singular  values.  We 
also  show  these  singular  values  as  function  of  nonlinearity  strength  ao  in  Fig.  7 
(lower  panel).  We  first  note  that  the  larger  singular  value  of  the  ideal  variance 
response  to  the  changes  in  dissipation  grows  at  an  extremely  high  rate  as  a o 
grows  from  ao  =  0  to  ao  =  0.20.  Then,  as  ao  increases  further  up  to  ao  =  0.35, 
the  larger  singular  value  changes  insignificantly.  On  the  other  hand,  the  larger 
singular  value  of  the  qG-FDT  variance  response  to  the  changes  in  dissipation 
has  a  less  rapid  but  steady  growth  as  ao  changes  from  ao  =  0  to  ao  =  0.35.  The 
smaller  singular  value  has  an  opposite  trend,  i.e.,  for  the  ideal  operator  it  grows 
with  a  smaller  rate  than  for  the  qG-FDT  operator. 

Finally,  we  consider  the  monthly  response  of  the  variance  to  changes  in 
external  forcing  across  all  months  for  the  resonant  forcing  case  with  a o  =  0.35 
(see  fourth  raw  in  Table  1).  This  is  an  extremely  difficult  test  case  for  the 
quasi-Gaussian  approximation  since  a  purely  Gaussian  approximation  would 
give  zero  response  [25] .  The  largest  singular  value  of  the  ideal  response  for  each 
month  is  reported  in  Fig.  8.  There  is  an  extremely  strong  variance  response  at 
some  months  and  a  weaker  one  elsewhere.  This  difficult  trend  is  captured  very 
well  by  the  qG-FDT  operator  (see  Fig.  8).  The  next  largest  singular  value  for 
all  months  is  almost  two  orders  of  magnitude  smaller  and  the  smallest  one  is 
identical  zero;  both  of  them  are  well  captured  by  qG-FDT  but  not  reported  here. 
This  provides  unambiguous  evidence  that  the  qG-FDT  algorithm  produces  the 
correct  correlations  among  the  modes  with  high  skill. 
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6.  Conclusions 


Time-periodic  FDT  is  tested  in  this  paper  on  a  stochastic  nonlinear  triad 
model  with  periodic  forcing  and  exactly  solvable  first  and  second  order  statistics. 
The  model  is  designed  to  mimic  the  seasonal  cycle  in  the  system  of  two  Rossby 
waves  nonlinearly  interacting  with  a  zonal  jet  where  the  Rossby  waves  are  forced 
by  moist  processes  and  the  zonal  jet  is  driven  by  polar  temperature  gradient. 

We  provided  a  brief  introduction  to  the  general  time-periodic  FDT,  in  Sec¬ 
tion  2.  Then,  in  Section  3,  we  discussed  the  stochastic  triad  model  and  showed 
how  to  find  analytically  first  and  second  (and,  in  principle,  any)  order  statistics. 
Next,  we  explained  how  these  exact  statistics  are  utilized  to  compute  the  ideal 
response  of  the  system  for  the  tinre-averaged  mean  and  variance  to  the  pertur¬ 
bations  of  external  forcing  or  dissipation.  The  ideal  response  is  then  used  as  a 
benchmark  for  the  response  operators  obtained  via  FDT.  Then,  in  Section  4, 
we  showed  how  to  construct  quasi-Gaussian  approximation  to  the  operators  for 
time-averaged  mean  and  variance  response  to  the  changes  in  forcing  and  dissi¬ 
pation  for  the  triad  model  (35),  (36).  Here,  we  also  need  to  utilize  the  exact 
mean  and  covariance  in  order  to  construct  the  Gaussian  pdf  with  the  same  first 
and  second  order  moments  as  in  the  original  system.  Then,  in  Section  5,  we 
develop  a  series  of  stringent  and  unambiguous  tests  for  the  qG-FDT  using  the 
triad  model  (35),  (36),  where  we  compute  the  skill  of  the  qG-FDT  operators  for 
the  mean  response  to  the  changes  in  the  forcing  and  for  the  variance  response 
to  the  changes  in  dissipation  with  the  corresponding  ideal  operators.  We  have 
found  that  in  the  near-Gaussian  regime,  qG-FDT  has  very  high  skill  for  the 
mean  response  to  the  changes  in  the  forcing.  Moreover,  surprisingly  this  skill 
stays  high  even  when  the  system  departs  from  Gaussian  regime  and  even  has 
two  peaked  pdf  as  in  Fig.  5.  On  the  other  hand,  qG-FDT  is  found  to  have  less 
skill  in  approximating  the  ideal  variance  response  to  the  changes  in  dissipation. 
We  attribute  this  behavior  of  qG-FDT  to  the  fact  that  the  variance  response 
to  the  changes  in  dissipation  is  computed  using  the  fourth  order  two-time  cor¬ 
relation  function  where  non-Gaussian  effects  can  be  more  significant  while  the 
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mean  response  to  the  changes  in  forcing  is  computed  using  the  second  order 
two-time  correlation  function.  Moreover,  we  considered  three  different  averag¬ 
ing  times,  i.e. ,  averaging  over  one  month,  one  season,  and  the  whole  year.  In 
our  test  cases,  the  pdfs  averaged  over  longer  time  (season  or  year)  were  closer 
to  Gaussian  state  than  the  pdfs  averaged  over  shorter  time  (month  or  season). 
This  lead  to  higher  skill  of  qG-FDT  for  the  statistics  averaged  over  the  whole 
year. 

There  are  a  number  of  directions  to  proceed  with  the  future  work.  One 
can  consider  a  more  general  type  of  perturbations  of  external  forcing,  SF  = 
a(u,  s)6f(t ),  where  a(u,  s)  actually  depends  on  s  as  in  the  case  of  the  perturba¬ 
tions  of  the  amplitude  of  external  forcing.  Another  interesting  question  to  ask 
is  whether  we  can  replace  the  response  function  R(t)  given  by  Eq.  (11)  with  the 
corresponding  correlation  only  within  the  averaging  time  interval.  This  may 
seem  a  simpler  approach  but  it  does  not  have  solid  theoretical  justification. 
One  of  the  natural  applications  of  the  time-periodic  FDT  is  to  assess  the  ef¬ 
fects  of  seasonal  cycle  on  climate  change  using  information  theory  for  finding 
the  most  “dangerous”  perturbations  [23,  27].  Also  to  make  an  approximation 
for  the  FDT  with  the  skill  beyond  qG-FDT,  one  can  use  moment  estimators 
as  discussed  in  [23,  26,  1].  Finally,  the  time-periodic  FDT  developed  in  [27] 
and  studied  here  can  be  applied  to  realistic  General  Circulation  Models  with 
seasonal  cycle  to  address  practical  climate  change  issues.  The  results  developed 
here  for  the  exactly  solvable  triad  model  should  provide  important  guidelines 
for  the  behavior  of  time-periodic  FDT  for  these  more  complex  systems.  In  par¬ 
ticular,  the  model  also  allows  for  exact  statistical  solutions  with  time-periodic 
dissipation  and  other  straightforward  more  general  perturbations. 
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A.  Second  and  third  order  statistics  and  ideal  variance  response  to 
the  change  in  forcing  and  dissipation 

We  show  how  to  find  the  second  order  statistics  of  u2  in  the  time-periodic 
equilibrium  state  of  the  triad  system  (35),  (36).  By  definition,  we  have 

Var(u2(t))  =  (\u2(t)\2)-\(u2(t))f.  (71) 


=  -72(2  t-s-r) 


(i/j(s,r))f2(s)fZ(r)dsdr 


We  use  Eq.  (42)  to  find 

<m*)i2>  =  e-2^-*o)(|U2oi2)+  r  f 

J to  J 

+e-72(t-to)  e-l2(t-s)  tty*  (s,t))f2(s)*ds  +  c.c. 

+  ^( 1-e-2-^ ). 


(72) 


Next,  we  consider  the  time-periodic  equilibrium  regime  by  taking  the  limit  to 


— oo 


{\u2(t)\'2)eq  =  +  f  [  e  72(24  s  r)(ip{s,r))f2(s)f2(r)dsdr. 
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Therefore,  the  equilibrium  variance  becomes 

Vareq(u2(t))  =  +  J  J  e-72(2*-s-r)^(S)?,)}  _  (ip(s,t))(tp{r,t))*^f2(s)f2(r)dsdr,  (73) 

where  ( ))  is  given  by  Eq.  (47).  Similarly,  we  find  the  cross-covariance  in 
time-periodic  equilibrium 


CoVeq(u2(t),U2(t)*)  =  ( U2(t )2}  -  ( U2(t ))5 
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where  [11,  12] 


—  g  i({J(s,t))+{J(r,t)))—^(Var(J(s,t))+Var(J(r,t))+2Cov(J(s,t),J(r,t))) 


and  ( J(s,t ))  and  Var(J(s,t))  are  given  by  Eq.  (48),  (49),  respectively,  and 

2  2 

Cov(J(s,  t),  J(r,  t))  =  —  fl  +  27i(max(s,  r)  —  t)  —  e~'1l('t~s'>  —  +  e_7lls_T,l  V  (76) 
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Finally,  we  find  the  cross-covariance  [11,  12] 

Coveq(u2,ui)  =  ^  3°  f  f2(s)e~72^~s)(ip(s,t))  (l  -  e~7l(t~s))ds.  (77) 

Similarly,  we  compute  the  triple  correlator  ^n\u\v!<^eq  for  u'j  =  Uj  —  {uj),  j  =  1,2 
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Below,  we  drop  the  subscript  “eq”  for  simplicity  of  notation  and  assume  that 
all  the  statistics  are  computed  in  the  time-periodic  equilibrium  regime.  Here  we 
show  how  to  find  the  mean  and  covariance  of  the  triad  (55)  in  the  time-periodic 
equilibrium  regime.  We  find  that 
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where  (iq)  is  given  in  Eq.  (39)  and  (u2)  is  given  in  Eq.  (46).  The  covariance 
matrix  of  the  triad  model  (35)  and  (36)  in  the  time-periodic  equilibrium  has  the 
form 
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where 
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Now,  we  obtain  the  ideal  variance  response  to  the  changes  in  forcing 
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where 
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and  Eqs.  (81)  were  also  used. 

Now,  we  obtain  the  ideal  variance  response  to  the  changes  in  dissipation 
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where 
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and  the  derivatives  are  computed  at  the  point  7  =  {h,i2)T  =  0.  We  use 
Eqs.  (73)  and  (74)  to  find 
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d(J(s,t )}  1  d  Tr  . 


=  -ao 


fds'  f  (s'  -s")h{s")e-^s'-s"Us", 

J  s  J — oo 

Var(J(s,t ))  +  — ^(i  -  s)(l  -  e_7l(t_s(V 
7i  7i  v  7 

/  d(J(s,t )}  d(J(r,t))  1  d  yr  , 

{ V’(M)V’M))  i - +i — — - -  —  Var(J(s,t)) 


d'yi 


d'yi  2  d'yi 


-|  O  A  ’ 

-^Q^Var(J(r,t))  +  —Cov(J(s,t),J(r,t)) 


Cov(J(s,t),  J(r,t))  =  — — Cov(J(s,t ),  J(r,t ))  —  aia?  (\  +  2(max(s,  r)  —  t) 

7i  27i  v 


d'yi 


Next,  we  compute 
dVar(u2 )  erf 

^72  2 7| 

and 

dCov(v,2(,  u 2) 

572 


+(t  -  s)e_7l(t_s)  +  (t  -  r)e"7l(t-r)  -  -  r|e-7l|s-r|)  (86) 


s,r))  ~  {i’(s,t))(ip(r,t))*)f2(s)fZ(r)dsdr,  (87) 


f  j  (2 t  -  s  -  r)e 

J  —  OO  J  —  OO 


[  [  (2  t-s-r)e  72(24  s  r)Uilj(s,t)ip(r,t))  -  (ip(s,t))(4>(r,t))'\  f2(s)f2(r)dsdr,  (88) 

J  —oo  J —00  '  ' 
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Case 

7i 

72 

O’  1 

02 

OJO 

a0 

/2O) 

1  Near-Gaussian 

0.1 

0.2 

0.1 

0.1 

107t/6 

1 

0.2ei(5,r/6t+1) 

2  Non- Gaussian 

season 

0.2 

0.4 

0.3 

0.1 

IOtt/6 

1 

gi(57r/6t+l)  y 

3  Non-resonant 

forcing 

0.1 

0.1 

0.2 

0.2 

7 r 

0.1,  0.2,  ...,0.5 

0.2ei(5,I'//6't+1) 

4  Resonant  forcing 

0.1 

0.1 

0.2 

0.2 

57t/6 

0.05,  0.10,  ...,  0.35 

0.2e*(5,I'/6't+1) 

Table  1:  Description  of  parameters  for  different  test-cases.  All  four  test-cases  have  the  same 
forcing  of  the  first  mode,  fi(t)  =  sin(57r/6  •  t). 


Table  2:  Ideal  and  qG-FDT  operators  for  the  time- averaged  mean  response  to  the  changes  in 
external  forcing  and  corresponding  singular  values  for  test-cases  1  and  2  from  Table  1. 
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0.23 

0.22 


Var(x, ) 

I  I  I  I . . . .  I 


0.21 


Cov(x2,x3) 


Var(x  ' 


1  2  3  4  5  6  7  8  9  10  11  12 


Figure  1:  (left  column  and  the  upper  two  panels  of  the  right  column)  First  and  second 
order  equilibrium  statistics  for  the  triad  (55)  with  the  parameters  given  by  the  second  row  in 
Table  1.  Analytically  obtained  statistics  (solid  line)  are  compared  with  the  results  of  Monte 
Carlo  averaging  (pluses),  (lower  two  panels  of  the  right  column)  Skewness  of  X2  computed 
through  Monte  Carlo  averaging  and  the  triple  correlator  (x\  x\  x'^)  computed  both  through 
Monte  Carlo  averaging  and  analytically.  Note  that  time  is  measured  in  months  here. 
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MONTH  1  SEASON  3  YEAR 


Figure  3:  Solid  line  shows  monthly,  seasonally,  and  annually  averaged  pdfs  for  the  near- 
Gaussian  test  case  with  non-Gaussian  season  given  in  the  second  row  of  Table  1;  the  pdfs  are 
rescaled  to  have  mean  zeros  and  variance  one.  Dashed  line  shows  standard  normal  distribution 
iV(0, 1)  (on  the  first  and  third  panels  the  solid  and  dashed  lines  are  on  top  of  each  other). 
The  pdfs  are  obtained  via  Monte  Carlo  simulation. 


ao 

Variance 

Skewness 

Flatness 

MONTH  1 

SEASON  3 

YEAR 

MONTH  1 

SEASON  3 

YEAR 

0.0 

0.100 

0.0020 

0.0015 

-0.0004 

3.0084 

3.0082 

3.0082 

0.1 

0.101 

-0.0015 

0.0009 

-0.0000 

3.0047 

3.0041 

3.0050 

0.2 

0.106 

-0.0009 

0.0068 

-0.0003 

3.0134 

3.0132 

3.0106 

0.3 

0.109 

-0.0239 

0.0498 

0.0001 

3.0192 

3.0467 

3.0494 

0.4 

0.129 

-0.0460 

0.1244 

0.0003 

3.1232 

3.1822 

3.2842 

0.5 

0.149 

0.0062 

0.1755 

0.0054 

3.5395 

3.4259 

3.6971 

Table  3:  Variance  of  X2  averaged  over  full  year  and  skewness  and  flatness  of  a? 2  averaged 
over  first  month,  third  season,  and  full  year  for  different  values  of  ao-  Non-resonant  forcing 
test-case  with  the  parameters  from  the  third  row  of  Table  1. 
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Figure  7:  (upper  panel)  Largest  (solid  line),  medium  (dashed  line),  and  smallest  (dotted  line) 
singular  values  of  the  mean  response  to  the  changes  in  external  forcing  as  functions  of  ao- 
Singular  values  of  the  ideal  operator  are  shown  with  circles  and  those  of  qG-FDT  operator 
are  shown  with  crosses.  The  response  operator  is  averaged  over  the  first  month.  Resonant 
forcing  test-case  with  the  parameters  from  the  fourth  row  of  Table  1  is  shown.  The  singular 
values  are  taken  from  Table  7.  (lower  panel)  Same  as  upper  panel  but  for  variance  response  to 
the  changes  in  dissipation.  The  larger  singular  value  is  shown  with  solid  lines  and  the  smaller 
singular  value  is  shown  with  dotted  lines.  The  singular  values  are  taken  from  Table  8 
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largest  singular  value  of  R, 


s 


figure  8:  Largest  singular  value  of  the  variance  response  to  the  chang> 
iach  month.  Thick  dashed  line  connects  the  singular  values  of  the  ide; 
;he  dotted  line  connects  the  corresponding  values  of  the  qG-FDT  resp< 
:orcinc  test-case  with  the  oarameters  from  the  fourth  row  of  Table  1 


ao 

Month 

Season 

Year 

Ideal 

FDT 

Ideal 

FDT 

Ideal 

FDT 

10.000 

9.787 

10.000 

9.960 

10.000 

9.919 

0.0 

0.317 

0.316 

0.317 

0.316 

0.317 

0.319 

0.317 

0.315 

0.317 

0.315 

0.317 

0.319 

10.025 

9.795 

10.020 

9.960 

10.000 

9.906 

0.1 

0.330 

0.328 

0.308 

0.307 

0.318 

0.320 

0.329 

0.328 

0.308 

0.305 

0.318 

0.320 

10.110 

10.062 

10.091 

10.212 

10.000 

10.084 

0.2 

0.343 

0.343 

0.299 

0.297 

0.319 

0.321 

0.340 

0.339 

0.297 

0.295 

0.319 

0.320 

10.285 

10.217 

10.234 

10.340 

10.000 

10.062 

0.3 

0.358 

0.356 

0.289 

0.286 

0.319 

0.321 

0.348 

0.346 

0.283 

0.280 

0.319 

0.321 

10.579 

10.568 

10.474 

10.639 

10.000 

10.112 

0.4 

0.373 

0.375 

0.277 

0.276 

0.320 

0.325 

0.352 

0.361 

0.265 

0.257 

0.320 

0.325 

10.957 

10.890 

10.782 

10.912 

10.001 

10.098 

0.5 

0.389 

0.397 

0.264 

0.257 

0.322 

0.330 

0.355 

0.377 

0.245 

0.229 

0.322 

0.329 

Table  4:  Singular  values  of  ideal  and  qG-FDT  operators  for  the  time-averaged  mean  response 
to  the  changes  in  external  forcing  for  the  third  test-case  from  Table  1  with  the  non-resonant 
forcing. 
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do 

Month 

Season 

Year 

Ideal 

FDT 

Ideal 

FDT 

Ideal 

FDT 

0.0 

2.020 

1.924 

2.020 

1.962 

2.020 

1.952 

1.428 

1.410 

1.428 

1.435 

1.428 

1.429 

0.1 

2.020 

1.976 

2.020 

2.008 

2.020 

2.002 

1.438 

1.439 

1.438 

1.461 

1.438 

1.457 

0.2 

2.022 

1.969 

2.022 

2.001 

2.022 

1.993 

1.475 

1.485 

1.474 

1.510 

1.474 

1.505 

0.3 

2.046 

2.118 

2.045 

2.153 

2.045 

2.144 

1.560 

1.549 

1.560 

1.578 

1.559 

1.570 

0.4 

2.252 

2.956 

2.242 

2.986 

2.242 

2.975 

1.653 

1.724 

1.650 

1.753 

1.647 

1.742 

0.5 

2.813 

4.958 

2.802 

4.998 

2.799 

4.977 

1.667 

1.827 

1.671 

1.870 

1.660 

1.850 

Table  5:  Singular  values  of  ideal  and  qG-FDT  operators  for  the  time-averaged  variance  re¬ 
sponse  to  the  changes  in  dissipation  for  the  third  test-case  from  Table  1  with  the  non-resonant 
forcing. 


a0 

Variance 

Skewness 

Flatness 

MONTH  1 

SEASON  3 

YEAR 

MONTH  1 

SEASON  3 

YEAR 

0.00 

0.10 

0.0076 

-0.0060 

0.0004 

3.0168 

3.0158 

3.0148 

0.05 

0.14 

-0.0901 

-0.0189 

-0.0000 

2.9984 

3.0033 

3.0256 

0.10 

0.25 

-0.2468 

-0.0947 

-0.0023 

2.7937 

2.8481 

3.1285 

0.15 

0.37 

-0.3052 

-0.1292 

-0.0047 

2.5346 

2.6302 

3.0800 

0.20 

0.47 

-0.3145 

-0.1309 

-0.0061 

2.3731 

2.4638 

2.8673 

0.25 

0.53 

-0.3019 

-0.1164 

-0.0072 

2.3008 

2.3626 

2.6852 

0.30 

0.57 

-0.2673 

-0.1017 

-0.0073 

2.2434 

2.2972 

2.5393 

0.35 

0.59 

-0.2310 

-0.0953 

-0.0069 

2.2335 

2.2791 

2.4535 

Table  6:  Variance  of  X2  averaged  over  full  year  and  skewness  and  flatness  of  #2  averaged  over 
first  month,  third  season,  and  full  year  for  different  values  of  ao.  Resonant  forcing  test-case 
with  the  parameters  from  the  fourth  row  of  Table  1. 
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a  o 

Month 

Season 

Year 

Ideal 

FDT 

Ideal 

FDT 

Ideal 

FDT 

10.000 

9.911 

10.000 

10.030 

10.000 

10.030 

0.00 

0.375 

0.348 

0.375 

0.405 

0.375 

0.405 

0.375 

0.342 

0.375 

0.401 

0.375 

0.401 

13.648 

13.472 

13.090 

13.102 

10.000 

9.986 

0.05 

0.387 

0.365 

0.355 

0.381 

0.377 

0.384 

0.284 

0.261 

0.271 

0.296 

0.377 

0.384 

18.631 

18.424 

17.434 

17.487 

10.001 

10.074 

0.10 

0.401 

0.380 

0.341 

0.360 

0.380 

0.384 

0.215 

0.206 

0.195 

0.206 

0.380 

0.384 

21.294 

20.672 

19.772 

19.487 

10.002 

10.067 

0.15 

0.413 

0.395 

0.329 

0.340 

0.381 

0.380 

0.194 

0.187 

0.167 

0.179 

0.381 

0.380 

21.987 

20.757 

20.368 

19.542 

10.004 

10.061 

0.20 

0.421 

0.417 

0.321 

0.328 

0.382 

0.382 

0.192 

0.193 

0.157 

0.170 

0.382 

0.382 

21.624 

19.738 

20.027 

18.608 

10.006 

9.969 

0.25 

0.428 

0.435 

0.315 

0.316 

0.382 

0.383 

0.198 

0.206 

0.157 

0.170 

0.382 

0.382 

20.815 

18.886 

19.294 

17.840 

10.008 

10.105 

0.30 

0.433 

0.447 

0.310 

0.307 

0.382 

0.382 

0.208 

0.226 

0.161 

0.175 

0.382 

0.381 

19.873 

17.880 

18.448 

16.946 

10.010 

10.177 

0.35 

0.437 

0.460 

0.307 

0.304 

0.383 

0.384 

0.220 

0.245 

0.167 

0.181 

0.382 

0.382 

Table  7:  Singular  values  of  ideal  and  qG-FDT  operators  for  the  time-averaged  mean  response 
to  the  changes  in  external  forcing  for  the  fourth  test-case  from  Table  1  with  the  resonant 
forcing. 
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ao 

Month 

Season 

Year 

Ideal 

FDT 

Ideal 

FDT 

Ideal 

FDT 

0.00 

2.020 

1.908 

2.020 

1.934 

2.020 

1.929 

1.428 

1.379 

1.428 

1.397 

1.428 

1.389 

0.05 

4.242 

2.646 

4.166 

2.636 

3.931 

2.532 

1.942 

2.197 

1.942 

2.230 

1.941 

2.206 

0.10 

9.917 

3.480 

9.726 

3.495 

9.159 

3.407 

1.939 

2.155 

1.939 

2.148 

1.936 

2.115 

0.15 

14.338 

4.228 

14.155 

4.287 

13.644 

4.264 

2.014 

2.252 

2.001 

2.237 

1.959 

2.098 

0.20 

16.643 

4.780 

16.514 

4.833 

16.171 

4.812 

2.126 

3.154 

2.090 

3.044 

1.983 

2.554 

0.25 

17.416 

5.654 

17.331 

5.693 

17.119 

5.546 

2.201 

3.225 

2.149 

3.077 

2.001 

2.387 

0.30 

17.317 

7.232 

17.263 

7.213 

17.132 

7.027 

2.227 

3.268 

2.169 

3.062 

2.012 

2.388 

0.35 

16.786 

8.437 

16.750 

8.485 

16.668 

8.305 

2.220 

3.103 

2.164 

2.959 

2.018 

2.392 

Table  8:  Singular  values  of  ideal  and  qG-FDT  operators  for  the  time- averaged  variance  re¬ 
sponse  to  the  changes  in  dissipation  for  the  fourth  test-case  from  Table  1  with  the  resonant 
forcing. 
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